{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates how to apply some filtering and processing on induction and gamma logs. The aim is to make the data more interpretable by removing features that are not related to the hydrogeology of the bore. For example, steel bolts in the bore casing give a large induction response but have no relation to the formation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import lasio # package for reading and writing las files\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import sys, os\n",
    "import datetime\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define key functions\n",
    "\n",
    "def plot_logs(ax, ser, title, logplot = True):\n",
    "    \"\"\"\n",
    "    @param ax: matplotlib axis object\n",
    "    @param ser: pandas series object containing values to be plotted\n",
    "    @param title: title of plot\n",
    "    @param logplot: boolean: if true the values are plotted in log10 space\n",
    "    \"\"\"\n",
    "    \n",
    "    ax.plot(ser.values, ser.index)\n",
    "    if logplot:\n",
    "        ax.set_xscale('log')\n",
    "       \n",
    "    ax.grid(which='major', linestyle='-', linewidth='0.5', color='grey')\n",
    "    ax.grid(which='minor', linestyle=':', linewidth='0.5', color='grey')\n",
    "    ax.set_title(title)\n",
    "\n",
    "\n",
    "def combination(n, r):\n",
    "    \"\"\"\n",
    "    Function for finding combinations of a sequene\n",
    "    @param: n. \n",
    "    \"\"\"\n",
    "    return int((math.factorial(n)) / ((math.factorial(r)) * math.factorial(n - r)))   \n",
    "\n",
    "def pascals_triangle(rows):\n",
    "    result = [] \n",
    "    for count in range(rows):\n",
    "        row = []\n",
    "        for element in range(count + 1): \n",
    "            row.append(combination(count, element))\n",
    "        result.append(row)\n",
    "        # count += 1 # avoidable\n",
    "    return result[rows-1]\n",
    "\n",
    "def binom_filter(x, kernel):\n",
    "    \"\"\"\n",
    "    Function that applies the binomial filter\n",
    "    @param x: 1D array on which to apply the filter\n",
    "    @param kernel: 1D array with binomial kernel\n",
    "    retrurn\n",
    "    \"\"\"\n",
    "    return np.mean(np.convolve(x, kernel, 'same'))\n",
    "    \n",
    "def run_filter(series, window, min_periods, filter):\n",
    "    \"\"\"\n",
    "    Function that applies the functoin to the pandas series. Mostly based on \n",
    "    https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.rolling.html\n",
    "    @param series: pandas series with the to be filtered data\n",
    "    @param: window: integer with filter window size\n",
    "    @param: min_periods:minimum observation to have a value returned\n",
    "    @param: filter: 'median', 'binomial' or 'mean': define what type of filter to use\n",
    "    \n",
    "    returns:\n",
    "    a Window or Rolling sub-classed for the particular operation\n",
    "    \"\"\"\n",
    "    if filter == 'median':\n",
    "        return series.rolling(window = window, min_periods = min_periods).median()\n",
    "    \n",
    "    elif filter == 'binomial':\n",
    "        kernel = pascals_triangle(window)/np.sum(pascals_triangle(window))\n",
    "        return series.rolling(window = window, min_periods = min_periods).apply(binom_filter, args = (kernel,), raw=True)\n",
    "    \n",
    "    elif filter == 'mean':\n",
    "        return series.rolling(window = window, min_periods = min_periods).mean()\n",
    "    \n",
    "    else:\n",
    "        print(\"Filter must be one of  'median', 'binomial' or 'mean'\")\n",
    "\n",
    "        \n",
    "        \n",
    "def filter_log(infile, bore_id, cols, windows, min_periods, csv_dir, fig_dir,\n",
    "               filters  = ['median'], remove_negs = True, log_space = True,\n",
    "               min_depth = np.nan, max_depth = np.nan):\n",
    "    \"\"\"\n",
    "    Function that filters a geophysical log based on the kwarg filter types and filtering\n",
    "    parameters.\n",
    "    \n",
    "    @param: infile: string: path to las file with geophysical log\n",
    "    @param: bore_id: string: borehole identification\n",
    "    @param: columns: comma delimited string with las file columns to filter\n",
    "    @param: window: integer or list/ array of integers with filter window size\n",
    "    @param: min_periods: integer or list/ array of minimum observation to have a value returned\n",
    "    @param: csv_dir: string: path to directory into which the csv with values will be written\n",
    "    @param: fig_dir: string: path to directory into which the filtered log figure will be written\n",
    "    @param: filters: sequence with 'median', 'binomial' or 'mean': define what type of filter to use\n",
    "    @param: remove_negs: boolean: if True then negative values will be removed prior to filtering\n",
    "    @param: log_space: boolean: if True then observations will be transformed into log space for filtering\n",
    "    @param: min_depth: float: if not nan then values shallower than this values are removed prior to filtering\n",
    "    @param: max_depth: float: if not nan then values deeper than this values are removed prior to filtering\n",
    "    \n",
    "    This function writes the filtered data into a csv, plots the filtered log and writes the filtering metadata into\n",
    "    a text file\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    # Check that windows and min period are list or array like\n",
    "    if not isinstance(windows, (list, tuple, np.ndarray)):\n",
    "        if isinstance(windows, int):\n",
    "            windows = [windows]\n",
    "        else:\n",
    "            print('Invalid dtype for windows. Please make an integer or list.')\n",
    "    \n",
    "    if not isinstance(min_periods, (list, tuple, np.ndarray)):\n",
    "        if isinstance(min_periods, int):\n",
    "            windows = [min_periods]\n",
    "        else:\n",
    "            print('Invalid dtype for min_periods. Please make an integer or list.')\n",
    "\n",
    "    if not isinstance(filters, (list, tuple, np.ndarray)):\n",
    "        if filters in ['median', 'binomial', 'mean']:\n",
    "            windows = [filters]\n",
    "        else:\n",
    "            print(\"Invalid entry for filter. Please make it string or list from ['median', 'binomial', 'mean']\")\n",
    "               \n",
    "    # Open the las file\n",
    "    las = lasio.read(infile)\n",
    "\n",
    "    df_logs = las.df()\n",
    "    \n",
    "    # GEt prefix for file naming\n",
    "    \n",
    "    prefix = infile.split('\\\\')[-1].replace('LAS','las').split('.las')[0]\n",
    " \n",
    "    # Iterate through columns\n",
    "    \n",
    "    columns = [x.strip() for x in cols.split(',')]\n",
    "    \n",
    "    for item in columns:\n",
    "        \n",
    "        # Remove any white spaces\n",
    "        item = item.strip()\n",
    "        \n",
    "        # Create a data series from the column of interest\n",
    "        series = df_logs[item]\n",
    "        \n",
    "        # Remove nans\n",
    "        \n",
    "        series = series.dropna()\n",
    "\n",
    "        # Remove negative values\n",
    "        if remove_negs:\n",
    "         \n",
    "            series = series[series>0]\n",
    "        # Now filter using parameters described above\n",
    "        \n",
    "        filtered = series.copy()\n",
    "        \n",
    "        for i, filter in enumerate(filters):\n",
    "        \n",
    "            filtered = run_filter(filtered, windows[i],\n",
    "                                  min_periods[i], filter)\n",
    "     \n",
    "        \n",
    "        # Clip the filtered if values are assigned\n",
    "\n",
    "        if np.isfinite(min_depth):\n",
    "            filtered =  filtered[filtered.index > min_depth]\n",
    "        \n",
    "        if np.isfinite(max_depth):\n",
    "            filtered =  filtered[filtered.index < max_depth]\n",
    "        \n",
    "        # Create a plot\n",
    "    \n",
    "        fig, (ax1,ax2) = plt.subplots(1,2, figsize = (8,8), sharex = True, sharey = True)\n",
    "        \n",
    "        # Plot the unfiltered data\n",
    "        plot_logs(ax1, series, logplot = logplot, title = \"unfiltered\")\n",
    "        ax1.invert_yaxis()\n",
    "        plot_logs(ax2, filtered, logplot = logplot, title = \"filtered\")\n",
    "        ax1.set_ylabel('depth (m)')\n",
    "        ax1.set_xlabel('apparent conductivity (mS/m)')\n",
    "        ax2.set_xlabel('apparent conductivity (mS/m)')\n",
    "        # Create directories if need be\n",
    "        \n",
    "        for dir in [csv_dir, fig_dir]:\n",
    "            new_dir = os.path.join(dir, bore_id)\n",
    "\n",
    "            if not os.path.exists(new_dir):\n",
    "                os.mkdir(new_dir)\n",
    "        \n",
    "        plt.savefig(os.path.join(fig_dir, bore_id, prefix + '_' + item + '.png'))\n",
    "        \n",
    "        filtered.to_csv(os.path.join(csv_dir, bore_id, prefix + '_' + item + '.csv'))\n",
    "        \n",
    "        plt.close()\n",
    "        \n",
    "    # Write out metadata file\n",
    "    filename = os.path.join(csv_dir, bore_id) + '\\\\filtering_metadata.txt'\n",
    "        \n",
    "    with open(filename, 'w') as f:\n",
    "        s= 'Downhole logs filtered on the ' + str(datetime.datetime.now()) + '\\n'\n",
    "        s+= 'The filtered variables are {} .\\n'.format(cols)\n",
    "        \n",
    "        filter_string = ', '.join(filters)\n",
    "        window_string = ', '.join([str(x) for x in windows])\n",
    "        min_periods_string = ', '.join([str(x) for x in min_periods])\n",
    "        s+= 'Filtering using {} filter(s) with window size {} and minimum period of {}\\n'.format(filter_string, window_string, min_periods_string)\n",
    "        if np.isfinite(min_depth):\n",
    "            s+= 'The log were clipped above {} m\\n'.format(str(min_depth))\n",
    "        if np.isfinite(max_depth):\n",
    "            s+= 'The log were clipped below {} m\\n'.format(str(max_depth))\n",
    "        f.write(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we will filter some logs. Bring a csv with the filtering parameters requires by the filter_log function into jupyter lab\n",
    "\n",
    "infile = r\"C:\\Users\\PCUser\\Desktop\\EK_data\\Boreholes\\gamma_induction\\EK_induction.csv\"\n",
    "\n",
    "df = pd.read_csv(infile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'C:\\\\Users\\\\PCUser\\\\Desktop\\\\EK_data\\\\EK_Induction_gamma\\\\2017_tfd_files_GA_ALT_system\\\\17BP01I\\\\17BP01I_induction_down_relogged.las'"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['Infile'].iloc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\PCUser\\Anaconda3\\envs\\hydrogeol_utils\\lib\\site-packages\\ipykernel_launcher.py:185: FutureWarning: The signature of `Series.to_csv` was aligned to that of `DataFrame.to_csv`, and argument 'header' will change its default value from False to True: please pass an explicit value to suppress this warning.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File not found. Check this path  C:\\Users\\PCUser\\Desktop\\EK_data\\EK_Induction_gamma\\2017_tfd_files_GA_ALT_system\\RN029663\\RN29663_induction_down.las\n",
      "File not found. Check this path  C:\\Users\\PCUser\\Desktop\\EK_data\\EK_Induction_gamma\\2017_tfd_files_GA_ALT_system\\RN029665\\RN029665_induction_down.las\n",
      "File not found. Check this path  C:\\Users\\PCUser\\Desktop\\EK_data\\EK_Induction_gamma\\2017_tfd_files_GA_ALT_system\\RN029653\\RN029653_induction_down.las\n",
      "File not found. Check this path  C:\\Users\\PCUser\\Desktop\\EK_data\\EK_Induction_gamma\\2017_tfd_files_GA_ALT_system\\RN029662\\RN029662_induction_down.las\n",
      "File not found. Check this path  C:\\Users\\PCUser\\Desktop\\EK_data\\EK_Induction_gamma\\2017_tfd_files_GA_ALT_system\\13BP01D\\13BP01D_induction_down.las\n",
      "File not found. Check this path  C:\\Users\\PCUser\\Desktop\\EK_data\\EK_Induction_gamma\\2017_tfd_files_GA_ALT_system\\17BP07D\\17BP07D_induction_down.las\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n",
      "Header section Parameter regexp=~P was not found.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File not found. Check this path  C:\\Users\\PCUser\\Desktop\\EK_data\\EK_Induction_gamma\\2017_tfd_files_GA_ALT_system\\RN007412\\RN007412_induction_down.las\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Header section Parameter regexp=~P was not found.\n"
     ]
    }
   ],
   "source": [
    "# Batch mode using parameters from the spreadsheet\n",
    "for index, row in df.iterrows():\n",
    "\n",
    "    infile = row['Infile']\n",
    "    bore_id = row['Bore_ID']\n",
    "    filter = row['filter']\n",
    "    filters = [row['filter_1'], row['filter_2']]\n",
    "    Columns = row['Columns']\n",
    "    windows = [row['window_filter_1'], row['window_filter_2']]\n",
    "    min_periods = [row['min_periods_filter_1'], row['min_periods_filter_2']]\n",
    "    csv_dir = row['csv_dir']\n",
    "    plot_dir = row['plot_dir']\n",
    "    remove_negs = row['remove_negs']\n",
    "    logplot = row['logplot']\n",
    "    min_depth = row['min_depth']\n",
    "    max_depth = row['max_depth']\n",
    "    \n",
    "    if filter == 1:\n",
    "        try:\n",
    "            filter_log(infile, bore_id, Columns, windows, min_periods,\n",
    "                       csv_dir, plot_dir, filters, remove_negs, logplot,\n",
    "                       min_depth = min_depth, max_depth = max_depth)\n",
    "        except FileNotFoundError:\n",
    "            print('File not found. Check this path ', infile)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C:\\Users\\PCUser\\Desktop\\EK_data\\EK_Induction_gamma\\2016_las_files_DGRT_system\\Y10A_Ord40\\Y10A_Ord40_conductivity_down.las\n",
      "Ord40\n",
      "['median', 'binomial']\n",
      "INDUCTION_CALIBRATED\n",
      "[30, 31]\n",
      "[30, 31]\n",
      "C:\\Users\\PCUser\\Desktop\\EK_data\\Boreholes\\gamma_induction\\induction_output\\csv_files\n",
      "C:\\Users\\PCUser\\Desktop\\EK_data\\Boreholes\\gamma_induction\\induction_output\\pngs\n",
      "True\n",
      "True\n",
      "1.97\n",
      "14.16\n"
     ]
    }
   ],
   "source": [
    "# To demonstrate the input parameters we print the last row of inputs\n",
    "\n",
    "print(infile)\n",
    "print(bore_id)\n",
    "print(filters)\n",
    "print(Columns)\n",
    "print(windows)\n",
    "print(min_periods)\n",
    "print(csv_dir)\n",
    "print(plot_dir)\n",
    "print(remove_negs)\n",
    "print(logplot)\n",
    "print(min_depth)\n",
    "print(max_depth)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "hydrogeol_utils",
   "language": "python",
   "name": "hydrogeol_utils"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
